We are IntechOpen, the world's leading publisher of Open Access books Built by scientists, for scientists

Open access books available 5,300

130,000 155M

International authors and editors

Downloads

Our authors are among the

most cited scientists 154 TOP 1%

Selection of our books indexed in the Book Citation Index in Web of Science™ Core Collection (BKCI)

# Interested in publishing with us? Contact book.department@intechopen.com

Numbers displayed above are based on latest data collected. For more information visit www.intechopen.com

# **Remote Sensing Studies of Urban Canopies: 3D Radiative Transfer Modeling**

Lucas Landier, Nicolas Lauret, Tiangang Yin, Ahmad Al Bitar, Jean‐Philippe Gastellu-Etchegorry, Christian Feigenwinter, Eberhard Parlow, Zina Mitraka and Nektarios Chrysoulakis

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/63886

#### **Abstract**

Need for better understanding and more accurate estimation of radiative fluxes in urban environments, specifically urban surface albedo and exitance, motivates development of new remote sensing and three‐dimensional (3D) radiative transfer (RT) modeling methods. The discrete anisotropic radiative transfer (DART) model, one of the most comprehensive physically based 3D models simulating Earth/atmosphere radiation interactions, was used in combination with satellite data (e.g., Landsat‐8 observa‐ tions) to better parameterize the radiative budget components of cities, such as Basel in Switzerland. After presenting DART and its recent RT modeling functions, we present a methodological concept for estimating urban fluxes using any satellite image data.

**Keywords:** DART, model, radiative transfer, albedo, thermal exitance, remote sens‐ ing, urban canopies

# **1. Introduction**

In today's world, concerns about sustainable energy production, related environmental issues, and their impact on urban agglomerations and their citizens are of high importance. In this context, the functioning of urban environments, including climate‐alternating interactions between atmosphere and human civilization, is studied globally. Remote sensing is a power‐ ful tool that is increasingly used in such studies due to improvements achieved in sensor

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

technology, as well as modeling of remote sensing measurements with respect to three‐ dimensional (3D) radiative and energy budget. One of these studies is the H2020 project URBANFLUXES (http://urbanfluxes.eu), which aims to develop methods estimating the anthropogenic heat flux (Q<sup>F</sup> ) of urban environments by employing remote sensing data [1]. In other words, its goal is to estimate the impact of human urban activities on the energy budget of city using satellite images. Important parts of the surface energy balance computation are the 3D radiative budget components specifically as urban surface albedo including thermal exitance. However, no remote sensing model able to simulate accurately spatial distribution of urban spectral albedo and exitance has been previously available. Three conditions must be fulfilled in order to achieve an acceptable solution of an urban radiative budget simulation:

– The model must consider explicitly the 3D architecture of urban environments and simulate radiance images and radiative budget of urban environment. Hence, apart from physical modeling considerations, the model must be able to work with urban databases, including spatial information on vegetation and digital elevation model.

– The model must work within any atmospheric conditions and possibly with air pollution of an urban environment. This requires to model radiative transfer of both the atmosphere above and the air among urban objects.

– An operational methodology must allow calibrating outputs of the remote sensing model in terms of 2D distribution of albedo and exitance (i.e., to produce image outputs). This calibra‐ tion is important, because one cannot expect to have access to the optical properties of all urban surface elements, which vary in space (e.g., tiles of roofs have different reflectance values depending on their age) and time (e.g., wet and dry roofs will exhibit different anisotropy of their reflectance).

**Figure 1.** DART calibration with a remote sensing image (Landsat‐8) for computation of the urban surface albedo over the city of Basel, Switzerland.

Here, we present a 3D radiative transfer model, DART, that fulfils these requirements and its recent improvements for studying urban and natural Earth landscapes with remote sensing acquisitions. We present also the approach that was recently designed and implemented to

assess the spatial distribution of DART input parameters: optical properties of surface elements (e.g., roofs, streets, vegetation). **Figure 1** summarizes this approach. It leads to DART simulated albedo and exitance maps that are calibrated with real‐time satellite acquisition. Ideally, these maps have a spatial resolution that is equal to that of satellite images that are used for the calibration.

# **2. The DART model**

# **2.1. Presentation**

DART computes radiation propagation in the three‐dimensional (3D) Earth/atmosphere system in the entire optical domain from the visible to the thermal infrared parts of the electromagnetic spectrum (EMS) [2–6]. The Earth surfaces and the atmosphere are simulated as a three‐dimensional (3D) medium (**Figure 2**). For any urban and natural landscapes, DART simulates the 3D radiative budget and acquisitions by satellite and airborne imaging radio‐ meters and LIDAR scanners aboard of space and airborne platforms. The DART model, developed in the CESBIO Laboratory (www.cesbio.ups‐tlse.fr/fr/dart.htm) since 1992, can work with any 3D experimental landscape configuration (atmosphere, terrain geomorpholo‐ gy, forest stands, agricultural crops, angular solar illumination of any day, Earth/atmosphere curvature, etc.) and instrument specifications (spatial and spectral resolutions, sensor viewing directions, platform altitude, etc.).

**Figure 2.** DART cell matrix of the Earth/atmosphere system. The atmosphere has three vertical levels: upper (i.e., just layers), mid (i.e., cells of any size), and lower atmosphere (i.e., same cell size as the land surface). Land surface ele‐ ments are simulated as the juxtaposition of facets and turbid cells.

DART has been successfully employed in various scientific applications, including develop‐ ment of inversion techniques for airborne and satellite reflectance images [7–9], simulation of airborne sensor images of vegetation and urban landscapes [10], design of satellite sensors (e.g., NASA DESDynl, CNES Pleiades, CNES LIDAR mission project [11]), among others. DART forward simulations of vegetation reflectance were successfully verified by real measurements [12] and also cross‐compared against a number of independently designed 3D

reflectance models (e.g., *FLIGHT* [13], *Sprint* [14], *Raytran* [15]) in the context of the RAdiation transfer Model Intercomparison (RAMI) experiment [16, 17].

DART creates and manages 3D landscapes independently from the RT modeling (e.g., visible and thermal infrared spectroradiometers, LIDAR, radiative budget). This multi‐sensor functionality allows users to simulate efficiently radiative transfer products of the same landscape as being captured by various sensors. Major scene elements are as follows: urban features, trees, grass and crop canopies, and water bodies. A DART simulated tree is made of a trunk, optionally with branches created with solid facets, and crown foliage that is simulated either as a set of facets or as a set of turbid cells, with specific vertical and horizontal distribu‐ tions of leaf volume density. Trees of several species with different geometric and optical properties can be located within the simulated scene of any user‐defined size randomly or based on exact coordinates. Urban objects (houses, roads, etc.) contain solid walls and roofs built from triangular facets. Finally, water bodies (rivers, lakes, etc.) are simulated as facets of appropriate optical properties. DART can use external libraries to import and to some extent also edit landscape elements, digital elevation models (DEMs), and digital surface models (DSM) produced by other software or measured in field (e.g., translation, homothetic and rotational transformations). Most importantly, the imported and DART‐created landscape objects can be combined into virtual Earth scenes of user‐defined complexity. This allows importation of whole cities from urban databases provided by city councils and urban planners. The optical properties of landscape elements and their geometry, as well as and properties of atmosphere, are specified and stored in adjacent SQL databases.

Atmospheric cells are used to simulate attenuation effects of satellite at‐sensor radiance and also to model influence of the atmosphere composition on radiative budget of Earth surfaces. The atmosphere can be treated just as an interface above the simulated Earth scene or as a light‐ propagating medium above and also within the simulated Earth scene, with cell sizes inversely proportional to the particle density. These cells are characterized by their gas and aerosol contents and spectral properties (i.e., phase functions, vertical profiles, extinction coefficients, spherical albedo). These quantities can be defined manually or obtained automatically from an atmospheric database. DART contains a database that stores the properties of major atmospheric gases and aerosol parameters for wavelengths between 0.3 and 50 µm. In addition, external databases can be imported, for instance, from the AErosol RObotic NETwork (AERONET; http://aeronet.gsfc.nasa.gov/) or from the European Centre for Medium‐Range Weather Forecasts (ECMWF; http://ecmwf.int/). Atmospheric RT modeling includes the Earth/ atmosphere radiative coupling (i.e., radiation that is emitted and/or scattered by the Earth and backscattered by the atmosphere towards the Earth). It can be simulated for any spectral band within the optical domain from the ultraviolet up to the thermal infrared part of the electro‐ magnetic spectrum. The Earth/atmosphere coupling was cross‐compared and successfully validated [18, 19] with simulations of the MODTRAN atmosphere RT model [20].

### **2.2. Recent improvements of DART**

Set of improvements had to be recently implemented in the DART model in order to provide the optimal products required by the URBANFLUXES project.

## *2.2.1. Work preparation*

The urban information used for the URBANFLUXES project was provided in the form of urban databases (for an example, see **Figure 7**). The two following adaptations had to be introduced to interface the databases with DART: (1) the format of 3D objects (i.e., triangular facets) stored in a common file format "\*. obj", and (2) the information about the vegetation, which was simulated as a turbid medium. Turbid vegetation was created from a set of characteristics, which for each tree in the urban scene provided geographical coordinates, physical dimen‐ sions, and optical properties. Those characteristics were obtained partially from external sources or partially from databases already exiting in DART (e.g., from database of optical properties for urban elements and vegetation).

### *2.2.2. Modeling of in situ camera acquisitions*

In the frame of the URBANFLUXES project, in situ cameras are used for acquiring images of urban canopy. The major objective is to better assess the properties (i.e., spectral reflectance/ emissivity and thermodynamic temperature) of the urban elements that are observed. The images acquired by these sensors can be very useful for validating the approach that we devised in order to retrieve maps of optical properties from satellite images. In this context, we introduced the modeling of these in situ sensors into DART. The major difficulty was linked to the fact that these sensors (e.g., fish‐eye cameras) have a wide field of view (FOV). Hence, different elements of a scene are not viewed along the same direction, which strongly impacts the radiometry and geometry of the acquired images. It is interesting to note that most remote sensing models neglect this fact: They consider that sensors have an infinitely small FOV. In short, DART can know simulate in situ cameras at any location in the Earth landscape (natural or urban), with any view direction, either upward or downward. Some examples of simulated images are shown here: downward looking sensor (**Figure 3**) that is on top of an urban district (Toulouse, France) and sensors with upward and oblique view directions (**Figure 4**).

**Figure 3.** DART simulation of a fish‐eye camera acquisition above an urban district of Toulouse (France). Left: Red‐ green‐blue spectral composite in natural colors; right: a thermal infrared image.

**Figure 4.** DART simulation of in situ cameras: (a) Trees with leaves simulated as facets viewed by a camera with an upward looking direction, (b) and (c) trees with leaves simulated as turbid material in an urban environment viewed by a camera with a horizontal view direction.

Work continues in order to generalize the simulation of in situ sensors in case of any Earth surface element: fluids, atmosphere, water, etc., for any Earth scene configuration (i.e., isolated scenes and infinitely repetitive scenes with/without continuity of local digital elevation model). In addition to help in understanding and calibrating remote sensing measurements, the simulation of in situ sensors can be very helpful for a number of applications: to determine the optimal location and view direction of in situ sensors and to assess local atmosphere and pollution impact on sensor acquisitions.

### *2.2.3. Atmosphere database*

In relation to the in situ cameras implementation, the DART atmosphere database and its management in the model were improved to offer higher flexibility of DART when dealing with different atmospheric conditions, especially in case of available in situ and/or satellite measurements of atmosphere. This atmosphere database was originally derived from simu‐ lations of the MODTRAN atmosphere radiative transfer model. Its use with DART was already validated with MODTRAN simulations. However, the simulation of aerosols was not as accurate as for gases. The DART atmosphere database was therefore completed using trans‐ mittance spectra for scattering and absorption mechanisms, derived from the atmospheric model MODTRAN. An interesting feature of this improvement is a new possibility to specify gas and aerosol amounts within the urban scene, independently of the atmosphere character‐ istics above the considered environment. This will be of a great help for assessment of local atmospheric properties and pollution impact using in situ sensor acquisitions. This improve‐ ment was accomplished by importing HITRAN [21] line‐by‐line cross section database (with specified temperature and pressure) for thermal infrared spectral domain, as well as the MPI‐ Mainz [22] cross section database for visible and near‐infrared spectral domains.

**Figures 5** and **6** show comparisons of DART and MODTRAN simulations in the visible and near‐infrared and the thermal infrared wavelengths, respectively. Both results demonstrate that the DART update and the introduction of new atmospheric database, combined with an improved radiative transfer modeling approach, brings DART simulations of atmosphere very close to MODTRAN 5.1 simulations, which is very encouraging especially for simulating ac‐ curately the in situ sensors.

**Figure 5.** DART (red) vs. MODTRAN 5.1 (blue) simulations in short wavelengths (UV, VIS, near infrared). Gas model: US standard, aerosol model: Rural, visibility: 23 km. (a) Sun irradiance, (b) bottom of atmosphere (BOA) radiance, (c) top of atmosphere (TOA) reflectance (*ρground* =0.5).

**Figure 6.** DART (red) vs. MODTRAN (blue) in the long wavelengths (thermal infrared). Gas model: Tropical, aerosol model: rural, visibility: 23 km. (a) Path radiance calculated at top of atmosphere (TOA) of scattered + emitted fluxes from atmosphere. (b) Direct transmitted radiance from Earth to TOA. (c) Total TOA radiance (*Tground* =299.15 K). (d) TOA brightness temperature (*Tground* =299.15 K).

#### *2.2.4. Decomposition of a sensor image into images per type of scene element*

A common difficulty for analyzing in situ sensor images (i.e., radiance images) is assessment of radiance and area proportion per type of surface material (e.g., wall, roof, atmosphere) inside an image pixel. Knowledge of the different radiance components is valuable for the "iterative calibration" that calibrates DART with remote sensing images as presented in Section 3.1.2. Hence, DART has been improved to facilitate in addition to the original sensor radiance image

*LDART,Δλ(xDART, yDART, Ω<sup>v</sup> )* the per‐pixel simulations of radiance *LDART,Δλ, n(xDART, yDART, Ω<sup>v</sup> )* and cross section *σ<sup>n</sup> (xDART, yDART, Ω<sup>v</sup> )* images of each type of scene element n in discreet directions along the sensor viewing direction *Ω<sup>v</sup> )*.

# **3. Computation of urban albedo: remote sensing calibration of DART**

Bottom of atmosphere irradiance and exitance are essential terms in the 3D radiative budget of urban and natural landscapes. In the short‐wavelength spectral domain, the ratio of exitance and irradiance is simply the albedo. DART can compute these terms for any time, date, and spectral band Δ*λ*, by using urban database inputs and actual atmosphere conditions that can be derived from satellite images or data from meteorological centers (e.g., ECMWF) or in situ measurements (e.g., AERONET network). The urban database, presented to DART in a \*.obj file format, must contain all buildings of the considered area, as well as the relevant information about urban vegetation in that area (i.e., geographic locations and physical dimensions of trees). It also requires distinction between roofs, walls, and streets, which means each facet in the database has to be registered in one of these groups. DART uses this registration to assign optical properties of the different scene elements. An example of a 3D urban model in nadir and an oblique view, which was built based on the urban database of the city of Basel, is shown in **Figure 7**.

The accuracy of DART simulated albedo *ADART* ,Δ*<sup>λ</sup>* and exitance *MDART* ,Δ*<sup>λ</sup>* corresponds with accuracy of the presented optical properties. **Figure 8**, which shows various views of DART simulated images of Basel in natural colors, illustrates importance of specific optical properties. Roofs and walls appear with different colors. These colors depend on the roof and wall optical properties. A major problem is that it is practically impossible to know the spatial distribution of these properties accurately beforehand, due to their spatial variability. In addition, they vary also in time. In consequence, a major objective of the URBANFLUXES project is to derive maps of optical properties per type of Earth surface element. The approach relies on the use of atmospherically corrected remote sensing data (in our example, Landsat‐8 images), in order

to compute the 3D radiative budget accurately. Two calibration methods were developed to compute products albedo more accurately: (1) a direct and (2) an iterative calibration.

**Figure 8.** DART simulated images of Basel, Switzerland: (a) pushbroom image with zoom, (b) airborne camera image, (c) satellite RGB image, (d) satellite TIR image.

### **3.1. DART calibration with satellite images**

D*irect calibration* is a straightforward method that calibrates DART without assessing the optical properties of the surface elements per pixel [23], while *iterative calibration* uses an iterative procedure in order to spatially derive the optical properties of the scene elements. In theory, this approach gives more accurate results, but is computationally more complex. Both methods are demonstrated on the example of DART simulation of the city of Basel (Switzer‐ land) calibrated using a Landsat‐8 multispectral image, with band specifications presented in **Table 1**, corrected of atmospheric effects by DLR using the ATCOR atmosphere correction model [24].


**Table 1.** Landsat‐8 band specifications.

It is important to note that the image of the 6th Landsat‐8 band was omitted, because of its spectral coincidence with the atmospheric absorption bands. However, this omission has no impact on presented methods as they are designed to work with any type of satellite data, with a number of spectral bands ≥2 and with any spatial resolution.

### *3.1.1. Direct method*

The direct method consists of the following five steps:

**Step 1.** Reflectance images corresponding to the satellite bands used for calibration are simulated. For these simulations, we set the optical properties to a realistic, but not exact value, as this piece of information is unavailable. The spatial resolution of images is equivalent to the size (*xDART* , *yDART* ) of the DART voxels (in this example to 2.5 m).

	- **◦** The corrected image corresponds to the so‐called direct sun illumination. Hence, bottom of atmosphere (BOA) irradiance is direct, without diffuse component. The reflectance is noted ρdd,BOA.
	- **◦** The corrected image corresponds to the actual BOA sun (direct) and atmosphere (diffuse) irradiance. The reflectance is then noted ρsd,BOA. It would be ρhd,BOA (hemispheric‐direct) if the atmosphere irradiance was isotropic.

Here, the atmospherically corrected satellite image corresponds to the case ρsd,BOA. The corresponding DART simulated term of interest is noted:

> r *DART DART DART S S BOA S atm sat sat* ,Δ , l( *x y E L t* , , , , , , W W W W ( ) ( ) )

where—*xDART* , *yDART* : coordinates of a given point in the DART simulated scene,

— Ω*<sup>S</sup>* , Ω*sat*: sun and satellite view direction, respectively,


**Step 2.** Obtained DART image *ρDART, Δλ(xDART, yDART, Ω<sup>S</sup> , ES, BOA(Ω<sup>S</sup> ), Latm(Ω), Ωsat, tsat)* was georeferenced and resampled to the spatial resolution of the Landsat‐8 image (*xsat*, *ysat*) to ensure exact correspondence in size and geolocation.

**Step 3.** A calibration factor *K*Δ*<sup>λ</sup>* (*xsat*, *ysat*, *t sat*) is computed per DART pixel (*xsat*, *ysat*). The objective of this factor is to mitigate the use of approximate optical properties. We use the following equation:

$$K\_{\Delta l} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, t\_{\text{sur}} \right) = \frac{\rho\_{\text{solular}} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, \boldsymbol{\Omega}\_{5}, \boldsymbol{\Omega}\_{\text{un}} \right)}{\rho\_{\text{solular}}, \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, \boldsymbol{\Omega}\_{5}, \boldsymbol{\Omega}\_{\text{un}}, t\_{\text{sur}} \right)}. \tag{1}$$
  $\text{Step 4. The map ofurban spectral albedo } A\_{\Delta l} \text{ is calculated per spectral band } \Delta l \text{ with:}$  
$$A\_{\Delta l} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, \boldsymbol{\Omega}\_{5}, E\_{S, \text{ROA}, \Delta l} \left( \boldsymbol{\Omega}\_{5} \right), L\_{\text{un}, \Delta l} \left( \boldsymbol{\Omega} \right), t\_{\text{sur}} \right) = K\_{\Delta l} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, t\_{\text{sur}} \right). \tag{2}$$
  $A\_{\text{DMTM}} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, \boldsymbol{\Omega}\_{5}, E\_{S, \text{ROA}, \Delta l} \left( \boldsymbol{\Omega}\_{5} \right), L\_{\text{un}, \Delta l} \left( \boldsymbol{\Omega} \right), t\_{\text{sur}} \right),$ 

where *ADART* ,Δ*<sup>λ</sup>* is the albedo computed by DART, for spectral interval Δ*λ*:

$$\begin{split} \boldsymbol{A}\_{\text{DARTA},\boldsymbol{\lambda}} \left( \mathbf{x}\_{\text{DART}}, \boldsymbol{\mathcal{y}}\_{\text{DART}}, \boldsymbol{\Omega}\_{\text{S}}, \boldsymbol{E}\_{\text{S}\ \boldsymbol{\mathcal{B}}\boldsymbol{\mathcal{A}}} \left( \boldsymbol{\Omega}\_{\text{S}} \right), \boldsymbol{L}\_{\text{atm}} \left( \boldsymbol{\Omega} \right), \boldsymbol{t}\_{\text{sat}} \right) &= \\ \frac{\boldsymbol{\rho}\_{\text{dh}} \boldsymbol{E}\_{\text{S},\boldsymbol{\text{B}}\boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{A}}} + \left[ \boldsymbol{\rho}\_{\text{dh}\,\boldsymbol{\text{A}}\boldsymbol{\lambda}} \left( \boldsymbol{\Omega} \right) \boldsymbol{L}\_{\text{atm}\,\boldsymbol{\text{A}}\boldsymbol{\lambda}} \left( \boldsymbol{\Omega} \right) \cos \left( \boldsymbol{\theta} \right) d\boldsymbol{\Omega}}{\boldsymbol{E}\_{\text{S},\boldsymbol{\text{B}}\boldsymbol{\text{A}}\boldsymbol{\text{A}}} + \left[ \boldsymbol{L}\_{\text{atm}\,\boldsymbol{\text{A}}\boldsymbol{\lambda}} \left( \boldsymbol{\Omega} \right) \cos \left( \boldsymbol{\theta} \right) d\boldsymbol{\Omega}}. \end{split} \tag{3}$$

In (3), *θ* is the zenith angle corresponding to the viewing direction Ω.

**Step 5.** Desired urban surface albedo *A* is calculated as the integral of all the spectral albe‐ dos *A*Δ*<sup>λ</sup>* (*xsat*, *ysat*, Ω*<sup>S</sup>* , *E<sup>S</sup>* ,*BOA*,Δ*<sup>λ</sup>* (Ω*<sup>S</sup>* ), *L atm*,Δ*<sup>λ</sup>* (Ω), *t sat*) across the entire spectrum, weighted by the BOA irradiance *E<sup>S</sup>* ,*BOA*,Δ*<sup>λ</sup>* (Ω*<sup>S</sup>* ).

**Figure 9.** Landsat‐8 atmospherically corrected image for an NIR band (864.6 nm) over the city of Basel, acquired on the 24th April, 2015.

Here, we consider an application of the method to the city of Basel. **Figure 9** shows the atmospherically corrected Landsat‐8 band 5 (**Table 1**) image.

The original DART image of the same spectral band without calibration, at the simulation resolution of 2.5 m, is shown in **Figure 10**. This image has been georeferenced using the spatial information retrieved from the Landsat‐8 satellite image.

**Figure 10.** DART reflectance image for an NIR band (864.6 nm), over the city of Basel, corresponding to a Landsat‐8 image, with erroneous (but realistic) optical properties used as input.

Image in **Figure 10** was resampled to the Landsat‐8 native spatial resolution using ILWIS image processing software. After computation of the calibration factor (step 3) and the albedo per spectral band, we integrated the final broadband albedo product as shown in **Figure 11**. In **Figure 11**, the color bar indicates the albedo values, and the axes simply indicate the pixels positions.

**Figure 11.** Urban surface albedo image computed using the direct calibration method for the city of Basel, with a Landsat‐8 multispectral image.

## *3.1.2. Iterative method*

Iterative calibration method computes the actual optical properties for each type of element in the scene per pixel of the satellite image or per group of pixels, which provides the desirable spatial variability in optical properties of present elements. If N types of elements are present in a single pixel, the reflectance of this pixel is formed by a number of different optical properties, that is, by N unknowns. Therefore, a set of equations is required to resolve the contribution of every surface element in a final reflectance value, which involves the image per type of scene elements described in Section 2.2. Taking into account a number of pixels, one obtains a system of equations, where a number of equations are equal or greater than the number of unknowns.

Consecutive steps of this approach are presented below, with N being the number of types of scene elements present in the studied urban area and parameter k being the order of the iteration initialized at k = 1:

**Step 1.** Regular grid is laid over the satellite image, with a mesh size equal to M M satellite pixels and with *M* > *N* . This ensures that each cell of the mesh contains a number of satel‐ lite pixels equal or larger than N.

**Step 2.** DART simulates the total radiance image *L DART* ,Δ*<sup>λ</sup>* (*xDART* , *yDART* , Ω*sat*) of the scene and also the radiance images per scene element n *L DART* ,Δ*λ*,*<sup>n</sup>* (*xDART* , *yDART* , Ω*sat*) in the satellite viewing direction Ω*sat*. Consequently, each surface element n viewed by a DART pixel d along the satellite viewing direction Ω*sat* has a reflectance value equal to *ρn*,*<sup>d</sup>* ,Δ*<sup>λ</sup> k* (*xDART* , *yDART* ).

It is important to note that reflectance and radiance are proportional terms, linked by the BOA sun irradiance:

( ) , ( ) ,Δ , , . , , , , . *DART DART DART sat DART DART DART sat DART BOA L x y x y E* l l l p r D D W W = (4)

The first iteration k = 1 assumes that the reflectance of an urban elements *ρn*,*<sup>d</sup>* ,Δ*<sup>λ</sup> k* (*xDART* , *yDART* ) is constant across the entire DART scene. Its plausible value is provided by a reasonable first guess, taken, for instance, from a recent airborne spectral data. Next step is to compute the reflectance value *ρn*,*<sup>d</sup>* ,Δ*<sup>λ</sup> <sup>k</sup>* +1 (*xDART* , *yDART* ) of the urban surface elements, which we will be used by DART in the follow‐up iteration k + 1. Objective of each iteration is to bring the DART reflectance closer to the reflectance of a real satellite image. In a DART pixel d, mean irradi‐ ance of an element type n is given by:

$$E\_{\text{DART},\lambda,\eta,d}\left(\mathbf{x}\_{\text{DART}},\mathbf{y}\_{\text{DART}}\right) = \frac{\pi \, L\_{\text{DART},\lambda,n}\left(\mathbf{x}\_{\text{DART}},\mathbf{y}\_{\text{DART}}\right)}{\rho\_{\pi}^{k}\left(\mathbf{x}\_{\text{DATT}},\mathbf{y}\_{\text{DATT}}\right)} \cdot \frac{\Delta\mathbf{x}\_{\text{DATT}}\,\Delta\mathbf{y}\_{\text{DATT}}\cos\theta\_{\text{sat}}}{\sigma\_{\pi,d}\left(\mathbf{x}\_{\text{DATT}},\mathbf{y}\_{\text{DATT}},\ \ \ \ \ \ \ \ \ \mathbf{S}\_{\text{sat}}\right)},\tag{5}$$

where *θsat* is the zenith angle of the satellite viewing direction and *σn*,*<sup>d</sup>* (*xDART* , *yDART* , Ω*sat*) is the cross section of the element of type n, which is viewed by the DART pixel d along the satellite viewing direction Ω*sat*. It is also an output of the DART model. The ratio *∆xDART* .*∆ yDART* .*cosθsat σn*,*<sup>d</sup>* (*xDART* , *yDART* , Ω*sat*) had to be introduced because the radiance of any DART pixel is simulat‐ ed per effective square meter of the pixel area *∆xDART* .*∆yDART* .*cosθsat*, while the radiance of the element n is expressed per effective square meter of the area of the surface element itself *σn*,*<sup>d</sup>* (*xDART* , *yDART* , Ω*sat*). Both radiances are equal, if only a single element is present in the considered pixel.

**Step 3.** DART radiance *L DART* ,Δ*<sup>λ</sup>* (Ω*sat*) and element *L DART* ,Δ*λ*,*<sup>n</sup>* (Ω*sat*) images are resampled in the satellite image spatial resolution (Δ*xsat*, Δ*ysat*). If selecting an appropriate spatial resolu‐ tion of a DART simulation, one can make the assumption that each satellite pixel m is composed of an integer number *D* 2 of DART pixels. Therefore, for any given satellite pixel m, where d is the index of the DART pixels in m (*d* ∈ 1 *D* 2 ) and n is the index of the element type (*n* ∈ 1 *N* ), one can define:

**•** Radiance as:

$$L\_{DART\lambda,\lambda,m} \left( \chi\_{\text{sat}}, \mathcal{Y}\_{\text{sat}}, \ \ \ \ \Omega\_{\text{sat}} \right) = \frac{\sum\_{d=1}^{D^2} L\_{DART\lambda,\lambda,d}}{D^2},\tag{6}$$

**•** Radiance of scene element n as:

( ) 2 ,Δ , , 1 ,Δ , , , Ω 2 , , *sat D DART n d d DART n m sat sat L L x y D* l l <sup>=</sup> = å (7)

**•** Irradiance of elements of type n as:

$$E\_{\text{DARTA},\lambda,\eta,d} \left( \mathbf{x}\_{\text{sar}}, \mathbf{y}\_{\text{sar}} \right) = \frac{\sum\_{d=1}^{D^2} E\_{\text{DARTA},\lambda,\eta,d} \left( \mathbf{x}\_{\text{DART}}, \mathbf{y}\_{\text{DART}} \right) \sigma\_{\eta,d} \left( \mathbf{x}\_{\text{DART}}, \mathbf{y}\_{\text{DART}} \right)}{\sum\_{d=1}^{D^2} \sigma\_{\eta,d} \left( \mathbf{x}\_{\text{DART}}, \mathbf{y}\_{\text{DART}} \right)},\tag{8}$$

and

**•** Cross section of element n viewed by pixel m as:

Remote Sensing Studies of Urban Canopies: 3D Radiative Transfer Modeling http://dx.doi.org/10.5772/63886 241

$$\sigma\_{n,m}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}}\right) = \sum\_{d=1}^{D^2} \sigma\_{n,d}\left(\mathbf{x}\_{DART},\mathbf{y}\_{DART}\right). \tag{9}$$

For a first approximation, Eq. (8) can be written similarly as Eq. (5):

$$\begin{array}{c} E\_{\text{DART},\text{A},\text{y},\text{an}} \left( \mathbf{x}\_{\text{au}}, \mathbf{y}\_{\text{au}} \right) = \frac{\pi \, L\_{\text{DART},\text{A},\text{y}} \left( \mathbf{x}\_{\text{au}}, \mathbf{y}\_{\text{au}}, \,\,\Omega\_{\text{au}} \right)}{\int \mathbf{p}\_{\text{n,m}}^{k} \left( \mathbf{x}\_{\text{au}}, \mathbf{y}\_{\text{au}} \right) d\mathbf{x} \,\,\Omega\_{\text{au}} \,\,\Omega\_{\text{au}} \,\,\Omega\_{\text{au}} \,\,\Omega\_{\text{au}}}. \end{array} \\ \begin{array}{c} \Delta\mathbf{x}\_{\text{au}} \,\Delta\mathbf{y}\_{\text{au}} \,\cos\theta\_{\text{au}} \\ \sigma\_{\text{u},\text{u}} \left( \mathbf{x}\_{\text{au}}, \mathbf{y}\_{\text{au}}, \,\Omega\_{\text{au}} \right) \\ \hline \end{array} \tag{10}$$

In this expression, *ρn*,*<sup>m</sup> k* (*xsat*, *ysat*) is the reflectance of the scene element of type n in the satellite pixel m at iteration k, as computed in the last iteration or equal to the initial value of k = 1. It assumes that the reflectance value of each DART pixel d corresponding to the satellite pixel m is constant.

**Step 4.** Resampled DART radiance is compared with the satellite radiance of all considered image pixels. If the difference between the two is acceptable, the procedure is terminated, and the desired urban albedo is the one computed by DART in the iteration of its simulation. Similarly to the direct calibration method, the albedo product is computed from all the spectral albedos. However, if the difference between the two radiances is according to user require‐ ments too large, the computation enters in a next step.

**Step 5.** The actual calibration of DART for each cell of the mesh defined in step 1 is conduct‐ ed for all *M* <sup>2</sup> satellite pixels and in each cell u of the mesh. Since *M* <sup>2</sup> > *N* , the number of pixels m analyzed in each cell u is larger than the number of different surface element types. Therefore, a deconvolution could be applied to retrieve the optical properties of each surface element present in the studied cell. Each cell u contains *M* <sup>2</sup> satellite pixels m, leading to a system of *M* <sup>2</sup> equations verifying if the DART image and the satellite image are equal:

$$\sum\_{\mathbf{n}=1}^{N} \underbrace{L\_{\text{DAT},\lambda,\eta,\mathfrak{n}}}\_{\mathbf{n}\_{\text{min}}} \left( \mathbf{x}\_{\text{sn}}, \mathbf{y}\_{\text{sn}}, \mathbf{\Omega}\_{\text{sn}} \right) = \overline{L}\_{\text{sn},\lambda,\mathfrak{n},\mathfrak{n}} \left( \mathbf{x}\_{\text{sn}}, \mathbf{y}\_{\text{sn}}, \mathbf{\Omega}\_{\text{sn}} \right), \forall m \in \left[ 1M^{2} \right]. \tag{11}$$

Obviously, the verification is negative if the two images differ. At iteration k, the DART radiance values *L DART* ,Δ*λ*,*n*,*<sup>m</sup>* (*xsat*, *ysat*, Ω*sat*) are computed with inaccurate approximated optical properties *ρn*,*<sup>u</sup> k* (*xsat*, *ysat*). If we define *L DART* ,Δ*λ*,*n*,*<sup>m</sup>* ' (*xsat*, *ysat*, Ω*sat*) as the DART radiance value computed from *ρn*,*<sup>u</sup> <sup>k</sup>* +1(*xsat*, *ysat*), and if we accept the fact that radiance values are proportional to reflectance values [Eq. (10)], then we can write:

$$\frac{L\_{\rm{DART}\lambda,\boldsymbol{\lambda},\boldsymbol{\eta},\boldsymbol{m}}\left(\boldsymbol{\chi}\_{\rm{sat}},\boldsymbol{\mathcal{y}}\_{\rm{sat}},\boldsymbol{\Omega}\_{\rm{sat}}\right)}{\boldsymbol{\rho}^{k}\_{\boldsymbol{n},\boldsymbol{\mu}}\left(\boldsymbol{\chi}\_{\rm{sat}},\boldsymbol{\mathcal{y}}\_{\rm{sat}}\right)} = \frac{L^{\prime}\_{\rm{DART}\lambda,\boldsymbol{\lambda},\boldsymbol{\eta},\boldsymbol{m}}\left(\boldsymbol{\chi}\_{\rm{sat}},\boldsymbol{\mathcal{y}}\_{\rm{sat}},\boldsymbol{\Omega}\_{\rm{sat}}\right)}{\boldsymbol{\rho}^{k+1}\_{\boldsymbol{n},\boldsymbol{\mu}}\left(\boldsymbol{\chi}\_{\rm{sat}},\boldsymbol{\mathcal{y}}\_{\rm{sat}}\right)}.\tag{12}$$

Consequently, we can rewrite the system of Eq. (11) to a system of *M* <sup>2</sup> equations, in which the unknowns are the expected reflectance values *ρn*,*<sup>u</sup> <sup>k</sup>* +1(*xsat*, *ysat*):

$$\sum\_{n=1}^{N} \frac{\rho\_{n,n}^{k+1}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}})}{\rho\_{n,n}^{k}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}})} \cdot L\_{\text{DMT}\lambda, \lambda, \eta, m}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}}) = L\_{\text{sat}\lambda, \eta, m}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}}), \forall m \in \left[1\,\text{M}^{2}\right]. \tag{13}$$

Resolving this system for each cell u of the mesh leads to the new desired reflectance values *ρn*,*<sup>u</sup> <sup>k</sup>* +1(*xsat*, *ysat*), which is used in the follow‐up iteration starting in step 2. This system will find no solutions, if one of the satellite pixels of the cell u contains N types of elements, but all other pixels contain just a single unique element of the same type. In this case, the solution needs a new adaptation of the mesh grid size.

The iterative method has been executed on the same case as the direct method, but so far we retrieved the new optical properties after only a single iteration. We show in **Figure 12** the preliminary result assessing the operational functioning of the method. On the left, we see the spectral Landsat‐8 image used for the calibration, for the 5th spectral band (864.6 nm). In the center is the corresponding initial DART reflectance image, resampled to the satellite resolu‐ tion. On the bottom is the new DART reflectance image, with updated optical properties for each element, on the new grid.

**Figure 12.** (Top left) Landsat‐8 image used for the calibration (NIR band). (Top right) DART image with initial optical properties (NIR band), resampled in the satellite resolution. (Bottom) DART image with updated optical properties (NIR band).

After one iteration, the computed mean relative error between the DART image and the Landsat‐8 image went from *rel* =0.3519 to *rel*,*new* =0.0998. It improved by a factor of 3.5. A clear advantage of this method compared to the direct one is the actual computation of new optical properties, which we are able to use in the simulation of other satellite images, whereas the direct calibration only works for a single image. A better accuracy is also expected after several iterations.

## **3.2. Albedo derived from climatological data**

Alternative methodology to derive urban albedo images has been adapted for a case when no satellite data of the area of interest are available. DART images are, in this case, simulated with optical properties of surface elements calibrated from the last available satellite image and with available atmospheric illumination conditions of direct sun irradiance *E<sup>S</sup>* ,*BOA*,Δ*<sup>λ</sup>* and diffuse sky irradiance *Eatm*,Δ*<sup>λ</sup>* . This information can be obtained from several sources including ECMWF (www.ecmwf.int), Meteo France (www.meteofrance.com), and other in situ mete‐ orological sensor networks.

In this alternative approach, DART calculates the so‐called white sky and black sky albedos.

The white sky albedo *ADART* ,*white sky*, Δ*<sup>λ</sup>* (*xsat*, *ysat*, *E<sup>S</sup>* ,*BOA*(Ω*<sup>S</sup>* ) =0, *L atm* = *Eatm π* , *t sat*) is computed for an irradiance coming solely from the atmosphere, that is, without any direct sunlight contri‐ bution. The black sky albedo *ADART, white sky, Δλ(xsat, ysat, ES, BOA(Ω<sup>S</sup> )=0, Latm=0, tsat)* is computed for the direct sunlight only. The desired final albedo product is then computed as a combination of both components:

$$A\_{\Delta\lambda} \left( \chi\_{\rm sat}, \chi\_{\rm sat}, \Omega\_S, E\_{s, \rm BOA} \left( \Omega\_S \right), E\_{\rm am}, t\_{\rm sat} \right) = E\_{s, \rm BOA} \left( \Omega\_S \right). \tag{14}$$

$$A\_{\rm DART, black \, sky, \lambda \lambda} + E\_{\rm am}, A\_{\rm DART, white \, sky, \lambda \lambda}$$

The black sky albedo must be pre‐computed for a set of sun directions such that the black sky albedo of an actual configuration (i.e., sun direction/date) can be derived by interpolation on pre‐computed black sky albedos.

### **3.3. Computation of thermal exitance images**

The calibration of the DART model in the thermal infrared domain, using satellite thermal infrared images, follows the same kind of approach. Indeed, the method will rely on spatial‐ ly setting up systems of equations over groups of pixels in the satellite image, to solve for the desired parameters. However, an additional difficulty comes from the fact that there are now 2 unknowns for each single measurement: the temperature *T* and the emissivity of the urban surface elements. Furthermore, the treatment of satellite pixels from thermal imagery only give information on the value of the product × *L* (*T* ), where *L* (*T* ) is Planck's law for temperature (*T* ). Hence, separating the variables by simply taking more pixels to create equations for the system is not adequate. The adopted strategy is to consider 2 thermal infrared images that correspond to 2 close spectral bands, with the assumption that emissivity is the same for these 2

spectral bands. This approach leads to the determination of the thermodynamic temperature and emissivity per surface element of type n. Then by considering the other spectral bands, if the satellite image has more than 2 thermal infrared bands, it leads to the determination of the emissivity of each type of surface element in those spectral bands. In the URBANFLUXES project, we will work with satellite images that have several thermal infrared bands, as shown in **Table 2**.


**Table 2.** Landsat‐8 thermal infrared band specifications.

# **4. Conclusions and future perspectives**

An innovative methodology using a 3D RT modeling and satellite reflectance data to derive urban surface albedo images has been designed and implemented. The method was tested using an atmospherically corrected Landsat‐8 multispectral image of the city of Basel. Although the preliminary results are encouraging, the newly presented methods have to be still properly validated. Different validation approaches are foreseen for testing their robust‐ ness:


A comprehensive comparison of the direct and iterative method technical performances is also required. It is especially important for future transformation of the methodology into a standardized operational approach.

Another important development lies in finalization of the thermal infrared spectral domain calibration. Combination of both calibrations in the visible and thermal part of the electro‐ magnetic spectrum would lead to a complete 3D radiative budget produced by DART at the same time as the other RT products.

To summarize, an original methodology for deriving urban surface albedo images from remote sensing data without a need of in situ measurements was implemented, although ground measurements could improve the final results. Both calibrations are computationally opera‐ tional, but the results for the iterative method are considered as preliminary, due to the limited implementation of only a single iteration. Iterative method is expected to be more accurate than for the simpler direct one. Nevertheless, the direct calibration might prove to be easier to use and more robust in case of inaccurate urban datasets, which depends on spatial resolu‐ tion and geometric co‐registration of DART and satellite imagery. Finally, both methods rely on the fact that DART is a reliable RT model simulating both satellite images and the radia‐ tive budget of cities with an acceptable accuracy. In a perspective of increasing remote sensing data availability, arising from new satellite systems such as Sentinel‐2, these approaches are expected to play an important role in future surveys of cities and their dynamic climate.

# **Author details**

Lucas Landier1\*, Nicolas Lauret<sup>1</sup> , Tiangang Yin1,2, Ahmad Al Bitar<sup>1</sup> , Jean‐Philippe Gastellu-Etchegorry<sup>1</sup> , Christian Feigenwinter<sup>3</sup> , Eberhard Parlow<sup>3</sup> , Zina Mitraka<sup>4</sup> and Nektarios Chrysoulakis<sup>4</sup>

\*Address all correspondence to: lucas.landier@cesbio.cnes.fr

1 CESBIO – UPS, CNES, CNRS, IRD, Toulouse, France

2 CENSAM, Singapore‐MIT Alliance for Research and Technology, Singapore, Singapore

3 UNIBAS – Basel University, Basel, Switzerland

4 FORTH – Hellas, Heraklion, Greece

# **References**

